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ABSTRACT 

We study the transition to runaway expansion of an initially stalled core-collapse supernova shock. The neu- 
trino luminosity, mass accretion rate, and neutrinospheric radius are all treated as free parameters. In spherical 
symmetry, this transition is mediated by a global non-adiabatic instability that develops on the advection time 
and reaches non-linear amplitude. Here we perform high-resolution, time-dependent hydrodynamic simula- 
tions of stalled supernova shocks with realistic microphysics to analyze this transition. We find that radial 
instability is a sufficient condition for runaway expansion if the neutrinospheric parameters do not vary with 
time and if heating by the accretion luminosity is neglected. For a given unstable mode, transition to runaway 
occurs when fluid in the gain region reaches positive specific energy. We find approximate instability crite- 
ria that accurately describe the behavior of the system over a wide region of parameter space. The threshold 
neutrino luminosities are in general different than the limiting value for a steady-state solution. We hypothe- 
size that multidimensional explosions arise from the excitation of unstable large-scale modes of the turbulent 
background flow, at threshold luminosities that are lower than in the laminar case. 
Subject headings: supernovae: general — hydrodynamics — shock waves — instabilities 



1. INTRODUCTION 

Following collapse of the core in a massive star, a hydrody- 
namic shock is launched when the central region reaches nu- 
clear density. Significant energy losses due to neutrino emis- 
sion and dissociation of infalling nuclei drain this sh ock of 
thermal energy and cause it to stall (e.g., Bethe 1990). For 
slowly-rotating progenitors, the re-absorption of a small frac- 
tion of neutrinos carrying away the gravitational binding en- 
ergy of the protoneutron star is thought to re-energ ize this 
shock and trigger an explosion dBefhe & Wilsonll9 85). How- 
ever, ab-initio radiation-hydrodynamic simulations in spheri- 
cal symmet ry fail to produce this out c ome for stars that form 
iron cores dLiebendorfer et al.l 1200 it iRampp & Jankal 120021; 
Thom pson et al.l 120031; ISumivoshi et al.l 12005b ! Success is 
only obtained for progenitors at the lowest en d of the mass 
range leading to core-col lapse (~ 8 - 10M©, iKitaura et al.1 
2006; iBurrows et all2007al) . 

Increasing the dimensionality of simulations improves con- 
ditions for explosion by enabling non-spherical hydrody- 
namic instabilities. In axisymmetry (2D), large scale shock 
oscillations combined with neutrino-driven convection lead 
to late-time - albeit somewhat marginal - explosion s in rep- 
resentative stellar progenitors ( Ma rek & Jank a 2009 and ref- 
erences therein). A recent three-dimensional (3D) hydro- 
dynamic study found that the extra spatial dimension could 
make conditions even more favorable, reducing the neutrino 
luminosity ne eded to start an explos ion by ~ 20% relative to 
axisymmetry dNordhaus et alj|20!ob . The way in which hy- 
drodynamic processes combine to cause this decrease is not 
well understood at present, however. In particular, a change 
in the methods and approximations employed can erase this 
dimensionality effect, suggesting that it may manifest only 
in some region of p hysical and/or numerical parameter space 
dHanke et al.11201 ll) . Given the marginality of 2D models, it 
is important to identify the conditions under which this addi- 
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tional reduction in neutrino luminosity occurs. 

The purpose of this paper and its companion is to system- 
atically examine the hydrodynamic processes responsible for 
the transition to explosion in a stalled core-collapse super- 
nova shock. As a first step, we study here the spherically 
symmetric case. Taking advantage of the simplicity of the 
flow geometry and borrowing tools from stellar pulsation the- 
ory, we develop a framework for identifying the relevant pro- 
cesses involved. The spherically symmetric case is also re- 
lated (in a time-averaged sense) to three-dimensional mod- 
els that experience small scale convecti on with no significant 
large scale shock deformati ons (e.g., iNordhaus et al.l |20ToT: 
Wongwat hanarat et al.ll2010l) . 

Spherically symmetric explosions, arising from boosted 
neutrino luminosities, involve a transition from a stalled con- 
figuration into runaway expansion on timescales longer than 
the dynamical time. In a realistic setting, the problem is fully 
time-dependent, and a complete analy sis must incl ude the mu- 
tual feedback of several e ffects (e.g., [Ja nka 2001). As a way 
to facilitate the analysis, Burrows & Goshy (1993) approxi- 
mated the problem as a steady-state accretion flow, and identi- 
fied a set of control parameters that determine it. They conjec- 
tured that the transition from accretion to explosion involved a 
global instability of the flow, with an associated critical stabil- 
ity surface in the space of control parameters. State-of-the-art 
radiation-hydrodynamic simulations show that multidimen- 
sional explosions take place at times when the background 
flow changes slowly (|Burrows et all l2007bt iMarek & Jankal 
2009; iSuwaetail 12010). making a steady-state treatment a 
reasonable starting point to study dimensionality effects. 

Global instabilities of the steady-state flow have been iden- 
tified in several linear stab i lity studies ([Yamasaki & Yamadal 
120051; iFoglizzo et al.ll2007t lYamasaki & Yamadall2007l) . Us- 
ing a realistic equation of state ( EOS) and weak interac- 
tions, lYamasaki & Yamadal d2007) found that as the neu- 
trino luminosity is increased, both oscillatory and non- 
oscillatory modes become unstable. The presence of these 
modes can be found in spherically-symmetric simulations 
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dating back to t he early days of t he delayed neutrino 
mechanism (e.g., IWilson et al.l Il986t IBurrows et al.l [T995b 
Buras et alj|2006bt lOhnishi et all 120061; iMurphv & Burrows! 
2008t IFernandez & T hompson! I2009at iNordhaus et all 120 lot 
Hanke et afll201 ll) . The instability mechanism behind these 
modes, their nonlinear development, and their connection to 
traditional diagnostics for explosion conditions are not com- 
pletely understoo d at present, howev er. 

In particular, Burrows & Goshy (1993) found that se- 
quences of steady-state models that satisfy a constraint on 
the neutrino optical depth are possible only up to a limiting 
neutrino luminosity, for a gi ven mass accretion rate. T his 
limiting value was found by Murphy & Burrows! (12008b to 
be close to the point where spherica l ly sy mmetric explosions 
occur, and by lYamasaki & Yamadal d2005l) to lie at the crit- 
ical stability curve in the neutrino luminosity vs. mass ac- 
cretion rate plane. These results provide support for the as- 
sumption that explosion is indeed obtained at the li miting 
steady-state luminosity (e.g., Peicha & Thom psonl201 ll) . The 
use of a realistic EOS in the linear analysis yields, how- 
ever, a flow that becomes unstable for luminosities ~ 30% 
lower than the limiting value, for a specific choice of param- 
eters l Yamasaki & Yamada 20Q3). IFernandez & Thom pson 
(2009a) found, in turn, that unstable modes develop into 
explosions when parametric microphysics and steady-state 
boundary conditions are used in time-dependent simulations. 

A more widely used diagnostic for explosion is the ratio of 
the advection to h eating time in the gain region dJanka & Keill 
Il998t I Thom pson 2000), which is known to have predictive 
powe r in the spherically symmetric case dThompson et alj 
2005). If an understanding of the hydrodynamic processes 
leading to explosion - in any number of dimensions - is to 
be obtained, then the connection between the limiting steady- 
state luminosity, stability in spherical symmetry, and proven 
explosion diagnostics needs to be clarified. 

We attempt to shed light on this issue here by study- 
ing the properties of unstable modes of the stalled super- 
nova accretion shock, their nonlinear development, and the 
their transition to a runaway solution. Our approach in- 
volves time-dependent hydrodynamic simulations, using a 
realistic EOS and weak interactions. The neutrino radia- 
tion field is treated parametrically, as is the fluid accret- 
ing onto the stalled shock. We focus on the evolution of 
the flow in between the neutrino sphere and the shock when 
the parameters describing this flow are systematically varied. 
Such a parametric approach has been employed previously to 
study properties o f exploding models and the effects of di- 
mensionality (e.g.,|Ta nka & M ullerll 19961 lOhnishi et"aT]l2006l: 
Mur phy & Burrowsll2008l) . 

This paper is organized as follows. Section 2 describes the 
physical assumptions and numerical methods employed. Sec- 
tion 3 revisits the limiting neutrino luminosity of the steady- 
state solution and its relation to the neutrino optical depth. 
Section 4 examines the instability affecting the post-shock 
flow in the linear and non-linear phase by borrowing analy- 
sis tools from stellar pulsation theory. By studying the ener- 
getics of the flow, the nature of the instability cycle, approx- 
imate stability criteria, and saturation mechanisms are iden- 
tified. Section 5 probes the validity of the instability criteria 
over a wide region of parameter space. The influence of nu- 
merical resolution, boundary conditions, and other parameters 
is also investigated. Finally, Section 6 contains a summary of 
our results and a discussion of the implications for the multi- 
dimensional case. Readers not interested in technical details 



can start in this section, where references to the relevant fig- 
ures and equations is made. 

2. METHODS 
2.1. Physical Model and Approximations 

We are interested in hydrodynamic instabilities that medi- 
ate the onset of explosion in a stalled core-collapse supernova 
shock, when the quantities that determine the global character 
of the accretion flow - mass accretion rate, neutrino luminos- 
ity, and protoneutron star radius - begin to evolve slowly rel- 
ative to the instability timescalefl This phase starts around 
100 — 200 ms after core bounce (e.g., Liebendorf er et all 
12001b . Prior to that, the system is still undergoing tran- 
sient ^eyohifcn^^Ej^ept fo^jJrograitorsat the lighter mass 
end dKitaura et a l. 2006; Buras et al. 2006a), simulations with 
neutrino transport obtain explosions for repr esentative pro- 
genitors only after this transient phase is over (IBurrows et al.l 
2007b; M arek & Jankdl2009t ISuwa et al.ll2010l) . 

The background ac cretion flow then evolves on timescales 
fhk„ >0j-l s (e g., lB^thl[T990l: lLiebendorfer"eTaT1 1200TI: 
Buras et al. 2006b|). In contrast, the dynamical time, advec- 
tion time, and thermal time of the flow in between the neutri- 
nosphere and the shock are approximately 

t s ^2M~ l ^ 2 r 3 7 /2 ms (1) 
f a dv~ 10r 7 v^ ms (2) 
t th ^20M l3 r^Q-l t2l ms, (3) 

where M\ 3 is the mass enclosed within the shock in units of 
1.3M©, r-i is the radial distance from the center of mass in 
units of 10 7 cm, v r $ the radial velocity in units of 10 9 cm s -1 , 
and 2net.2i the net specific neutrino energy source term, in 
units of 10 21 erg g" 1 s~ l (~ 1 GeV per baryon). These 
timescales are usually shorter than fbke by a factor of at least 
several, making a steady-state model a reasonable background 
state to study global hydrodynamic instabili ties and the tran- 
sition to explosion (IBurrows & Gosh y 1993). 

In this study we adopt a parametric approach in that we do 
not compute the neutrino radiation field and boundary condi- 
tions starting from a stellar progenitor. Instead, we regard 
the quantities determining those processes as free parame- 
ters, and study the behavior of the system as these param- 
eters are varied. This approach has previously been used 
by a number of authors to study the hy drodynamic response 
of the stalled accretion shock (e.g., Blondin et al.l 120031; 
Ohnishi et al. 2006; Murphy & Burrows 2008; Iwakami et al. 
20081; IFernandez & Thompson! l2009aFlNordhaus et al.l l2010t 
Hanke et al.ll2011l) . 

In our model, we ignore the interior of the protoneutron 
star, as we are interested in the dynamics of the fluid be- 
tween the neutrinosphere and the shock. Instead, we estab- 
lish an spherical inner boundary at the radius of the forming 
neutron star. A time-independent neutrino flux is imposed at 
this surface. We consider only electron-type neutrinos and 
antineutrinos, as these are the only species that exchan ge en- 
ergy w ith matter with a significant cross-section (e.g. jBethel 
1990). For simplicity, we establish a single neutrinosphere 
at a time-independent radius and assume that neutrinos 
stream isotropically from this surface. In a more realistic 

2 The mass of the protoneutron star is another fundamental parameter of 
the flow, but it must vary slowly in neutron-star-forming supernovae. 
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treatment, the difference between the energy averaged elec- 
tron neutri no- and antineutrin osphere radius is of the order of 
10% (e.g jBuras etai]200 6b). 

Both neutrino species are assumed to have a Fermi-Dirac 
spectrum with zero chemical potential. To remain close to re- 
sults from more detailed radiation-hydrodynamic simulations, 
we set the neutrino luminosities of both species to be equal, 
but all ow them to have different neutrinospheric temperatures 
(e.g.,[janka 2001). This is achieved by allowing the normal- 
ization of the neutrino luminosity to vary freely relative to the 
neutrinospheric temperature and radius (AppendixlAl. 

The contribution to neutrino heating from the accretion lu- 
minosity is not included explicitly. To first approximation, the 
amount of heating required to start an explosion depends only 
on the total flux incident on the gain region, which lies above 
the cooling region, so both contributions can be absorbed by 
the core luminosity in steady-state. However, instability of 
the post-shock flow generates changes in the mass accretion 
rate, hence a non-trivial feedback from the accretion luminos- 
ity can alter the stability prope rties . We discuss the implica- 
tions of this approximation in $4.51 

The mass accretion rate is taken to be constant in time. Only 
the point-mass gravity of the protoneutron star is considered, 
as the self -gravity of the envelope is a < 10% correction to the 
dynamics. This gravitating mass remains constant in time. 

Starting from a steady-state accretion flow defined as above, 
we evolve the system for different values of the control param- 
eters that set the background solution. In particular, we study 
the evolution of sequences of models with different neutrino 
luminosities, for fixed values of the mass accretion rate and 
neutrinospheric radius. Our goal is to identify the hydrody- 
namic processes by which a quasi-steady state configuration 
transitions into an exploding solution, and to determine the 
conditions (quantified by the control parameters) for which 
this transition occurs. 

2.2. Equations and Coordinate System 

Spherical symmetry is assumed throughout this paper, with 
the origin at the center of the protoneutron star. The time- 
dependent system is described by the equations of mass, mo- 
mentum, energy, and lepton number conservation in spherical 
polar coordinates, with source terms due to the gravity of a 
point mass M and charged-current weak interactions: 

% + \^pv r ) = (4) 
at r l or 

dv r dv r 1 dp GM 

— +v r — + -^f + ^-=0 (5) 

at or p or r l 

d ^ 1 d < r in GM c* 

-^— + ^Tr( v r[pe+p]) + pv r — r =J? ne t (6) 
Ot r l or r L 

-^-+v r — = r net . (7) 

ot or 

We denote respectively by p, vy, p, G, e, and Y e , the mass 
density, radial velocity, total pressure, gravitational constant, 
specific energy, and electron fraction. Source terms due to 
weak interactions are calculated in tabular form, with details 
provided in Appendix [A] The relevant contributions are the 
net rate of heating minus cooling per unit volume „Sf net , and 
the net rate per baryon of electron minus positron generation 

Tnet- 

The system of equations d4]i-(|7]i is closed with an EOS that 
yields the relation p(p,ei nt ,Y e ), with e; nt = e- v 2 r /2 the specific 



internal energy. We use the mo del oflShen etalj (l 998 ), as im- 
plemented bv lO'Connor&Ottl (I2010IFL The high-density part 
of EOS does not make a significant difference at the densities 
considered here, p < 10 11 g cm" 3 . The important components 
are nucleons, alpha particles, and heavy nuclei in nuclear sta- 
tistical equilibrium, supplemented by photons and electron- 
positron pairs with an arbitrary degree of degeneracy and rel- 
ativity. 

2.3. Initial Conditions 

The initial condition for time-dependent simulations con- 
sists of a steady-state spherical accretion flow with a shock 
at some radial distance R s from the center, obtained by 
solving equations (|4]i-(|7| without time derivatives. So- 
lutions upstream and downstream of the shock are con- 
nected through the Rank ine-Hugoniot jump conditions (e.g., 
lLandau & Lifshitzlll987l) . 

Conditions upstream of the shock are determined by re- 
quiring that the velocity, entropy, and electron fraction of 
the supersonic flow have specified values at a fixed radial 
coordinate r v f. In order to obtain values that are in close 
agree ment with more realistic one-dimensional simulations 
(e.g., lLiebend6rfer"eTaTI l200lt iBuras et alJ 1200 6b). we set 
r v f = 100 km for all simulations. At this radius, we set the 
radial velocity to the local free-fall velocity, the entropy per 
baryon to 5ks, and the electron fraction to 1/2. These param- 
eters together with the mass accretion rate M, the gravitating 
mass M, and the shock radius R s yield the mass density, in- 
ternal energy, electron fraction, and velocity upstream of the 
shock, by integrating the steady-state equations from r v f to R s . 
The sensitivity of our results to these choices is examined in 



For a given set of parameters, the shock radius R s is 
obtained by iteratively solving for the downstream flow 
from a trial shock position to R„, until an additional con- 
straint is satisfied. Following previo us studies of steady-state 
core-collapse superno va flo ws (e.g.. iBurrows & Goshy||1993l: 
Yamasaki & Yamada l2(M 120061 120071: iPejcha & Thom pson 
201 ll) . this additional closure relation is obtained by requiring 
that the neutrino optical depth from R v to the shock 



R, 



K,,dr 



(8) 



is equal to 2/3. By default, we set k v to the effec- 
ti ve absorption coefficient of electron-type neutrinos K e s = 
\/ K abs(Ksc + ftab.s) (Janka 2001), where 



^abs 



1.96xl0 _7 r^ 4 pi 7 n cm 



(9) 



corresponds to charged-current absorption and 



k sc ~ 0.51 x lO- 7 T^ A p w (Y n + Y p ) cm" 1 (10) 

to elastic scattering with nucleons, to lowest order in the 
neutron-proto n mass difference over the neutrino energy (e.g., 
lBruennl!985l) . In equations (f9b-(fT0T>. T v 4 is the electron neu- 
trinospheric temperature in units of 4 MeV, pio the mass den- 
sity in units of 10 10 g cm" 3 , and {¥^,,1^} the number of neu- 
trons and protons per baryon, respectively. To test the sensi- 
tivity of our results to this prescription for the optical depth, 
we also compute solutions with pure absorption («„ = re a bs), 

3 Available athttp://stellarcollapse.org 
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FIG. 1. — Initial profiles of density (a), entropy (b), electron fraction (c), 
and net neutrino energy generation rate Jzfnet / p (d) as a function of radius, 
obt aine d by solving the steady-state version of equations fiJ-jT) as described 
in jj2.3l Parameters correspond to our fiducial sequence ( 32.51 . Curves shown 
correspond to neutrino luminosities at the approximate instability thresholds 
for oscillatory modes (blue) and non-oscillatory modes (red), plus the limit- 
ing luminosity for a steady-state solution (black). 

or total absorption plus scatter ing (k v = k r \, s + k sc ), thereby 
bracketing our default choice ( jj5.21 i. 

Figure Q] shows characteristic radial profiles of density, 
entropy, electron fraction, and net specific neutrino energy 
generation rate for a few neutrino luminosities, given M = 
O.3M s" 1 andR v = 30 km. 

2.4. Time-dependent Implementation 

We use FLASH3.2 dDubev et al.1 120091) to evolve the sys- 
tem of equations ©-10. The public version of the code 
has been modified to include the EOS implementation of 
IP' Connor & Ottl ( 120101) . weak interaction rates, and a grid of 
variable spacing. Details about the latter are provided in Ap- 
pendix [B] 

We locate the inner simulation boundary at the neutri- 
nospheric radius R„. The outer boundary r max is chosen to 
be 1000 km for most models, corresponding to approximately 
four times the radius where the nuclear binding energy of al- 
pha particles equals their gravitational binding energy, 



where M enc is the mass enclosed at a radius r a . Above this 
radius , the shock accelerates signific antly durin g an explosion 
(e.g., Fernandez & Thompson 2009a and ^4.41 i. 

We use two types of boundary condition at R v . Our default 
implementation fixes all variables in the ghost cells to their 
initial values, obtained by continuing the steady-state solution 
inside R„. This prescription fixes the mass-, momentum-, and 
energy fluxes leaving the domain. A truly reflecting bound- 
ary condition would cause the shock to drift outwards with 
time due to the accumulation of mass around R v , because the 
latter does not recede (as it would do with a more realistic 
treatme nt, e .g.. IScheck et al.ll2006l). [Fernandez & Thompson! 
( 2009b) and Fernandez & Thompson ( 2009a) were able to use 
a reflecting boundary condition with a fixed neutrinospheric 
radius because their parametric cooling function was very 
centrally concentrated, resulting in accumulation of mass in 
a few cells outside R v , and thus eliminating the shock drift 
effect. The cooling function used here has, in contrast, a shal- 
lower radial dependence. 

To test the influence of the default inner boundary condi- 
tion on our results, we also run a few models that allow an 
arbitrary amount of mass to leave the domain, while still pro- 
viding pressure support to the accretion flow. The density and 
radial velocity in the ghost cells inside R v are set to 

p(r) = p (r) + p(r l )-p (r l ) (12) 
v r (r)=(^) 2 y r (n), (13) 

where r is the radial position of a given ghost cell, r\ corre- 
sponds to the center of the first active cell outside the inner 
boundary, and the subscript zero labels the initial steady-state 
solution. All other variables are set to have zero gradient and 
are thus copied from the innermost active cell into the ghost 
cells. The additional factor multiplying the velocity accou nts 
for the radial convergence of the flow (Ohnishi et al. 2006). 

The radial cell sp acing Ar is ratioed (as in, e.g., 
IStone & Normanlll992l) : Ar M /An = C, with i denoting cell 
number and £ some positive real number greater than one. We 
space our cells logarithmically, with £ = (r max /R„) N ', where 
N r is the total number of cells, and the size of the cell closest 
to the inner boundary Ar^ =R„(£- 1) (AppendixlBt. Guided 
by convergence tests, we choose either 880 or 1760 cells for 
our runs, altho ugh s ome models are evolved at resolutions up 
to N r = 3200 ( 315.21 . The corresponding values of (£— 1) are 
0.4%, 0.2%, and 0.1%, respectively, for r max = 1000 km and 
R v = 30 km. 

To avoid problems with the Riemann solver whenever the 
interna l energy becomes neg ative, we implement the prescrip- 
tion of iBuras et al.l d2006bl) adding the nuclear binding en- 
ergy of alpha particles, heavy nuclei, and the rest mass energy 
of electrons to the internal energy during the hydrodynamic 
step. The mass fractions of alpha particles and heavy nuclei 
are then advected as passive scalars in between calls to the 
equation of state. In some cases where the Mach number of 
the upstream flow is > 5, we add an additional zero point of 
7 x 10 17 erg g" 1 to prevent the internal energy to be negative at 
large radii. This constant shift makes a negligible difference 
in the evolution of the shock. 

To prevent the system from reaching the lower limit on 
Y e allowed by the equation of state implementation, and to 
account for the decrease in neutrino flux around the neutri- 
nosphere in a very crude manner, we multiply all the neutrino 
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TABLE 1 
Model Sequences Evolved 





M(M e s" 1 ) 


L Ve (10 52 ergs" 1 ) 




20 km 


0.1-0.5 


2.0-14.0 


1760 


30 km 


0.1-1 


1.5 - 13.7 


880 


40 km 


0.1-1 


1.0-8.9 


880 


r 1/2 


0.1-1 


1.0-5.5 


880 


30 km 


0.3 


3.8-4.6 


240-3200 



source terms by a suppression factor 

f™ v (p) = e- (plpo) (14) 

where po = 10 11 g cm" 3 is a fiducial density, chosen so as to 
match the characteristic density at the neutrinosphere. Given 
that the effective neutrino optical depth depends chiefly on 
density (eq. (8)), the factor in equation ( fl4l i amounts to a re- 
duction ~ e~ r " (e.g., iMurphv & Burrowsll2008l) . The choice 
of po has some mild influence on the threshold for instabil- 
ity, but it doesn't sign ificantly affect the relative differences 
between models ( jj5.21 >. 

A small initial transient is produced in the form of an out- 
going sound wave from the inner boundary and downgoing 
entropy and sound waves from the shock. The former acts as 
an initial perturbation. 

2.5. Models Evolved 

We evolve sequences of models with varying neutrino lu- 
minosity L Ve (= Ljy e ), for a given mass accretion rate M and 
prescription for the neutrino spheric radius R u . Our fiducial 
sequence corresponds to M = 0.3M Q s" 1 and R v = 30 km, 
with other parameters fixed to values as described in 32.31 
This default sequence closely resembles the flow obtained 
at late times in the evolution of a \5Mp, pr ogenitor (e.g., 
iLiebendorfer et al.ll200ltlMarek & Janka1l2009l) . 

We explore a wider region of parameter space with a grid of 
constant R u and M sequences such that R v = {20,30,40} km 

and M= {0.1,0.3,0.5,0.75, 1}M Q s _1 , with our fiducial se- 
quence as one of its members. W e add another sequence, 
following iBurrows & Goshvl d 19931) . that relates the neutri- 
nospheric radius to the neutrino luminosity through the Fermi- 
Dirac distribution at constant neutr inos pheric temperature and 
zero chemical potential (equation IIA3I with N Ur = 1), namely 

1/2 

R u oc LJ, ■ In this case we use the same grid in M as in the 
constant R v models. All sequences are summarized in Ta- 
ble Q] 

3. ON THE LIMITING NEUTRINO LUMINOSITY OF THE 
STEADY-STATE SOLUTION 

IBurrows & Goshvl i 19931) found that steady-state solutions 
to the accretion shock problem in the supernova context 
can exist only up to a maximum value of the neutrino lu- 
minosity (the Burrows-Go shy limit hereafter). Recently, 
Peicha & Thompson (2011) have related this limit to a crit- 
ical value of the ratio of sound speed to free-fall speed in the 
postshock region, in analogy with isothermal accretion flows. 
Here we provide an alternative (but equivalent) explanation. 

The existence of a limiting neutrino luminosity for the 
steady-state solution is related to the requirement that the neu- 
trino optical depth have a fixed value (eq. (8)). This closure 
relation is used for self-consistency with the assumed neutrino 
radiation field, but has no significance from the purely hydro- 
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Qj , , , , 1 , 1 , 1 

50 100 300 500 

shock radius [km] 

FIG. 2. — Neutrino optical depth r„ (eq. (8)) as a function of shock ra- 
d ius R s , for parameters corresponding to our fiducial sequence of models 
( j|2.5t . Each curve is obtained by fixing the neutrino luminosity and in- 
tegrating the time-independent version of equations (4j-{7) from an arbi- 
trary shock position to R v , keeping the upstream flow fixed. Curves with 
decreasing amplitude correspond to electron neutrino luminosities L„ t = 
{3.5,4.05,4.35,4.55,4.74,5} X 10 52 erg s -1 , respectively. The thick blue 
and red curves correspond to the approximate ins tabil ity threshold for oscil- 
latory and non-oscillatory modes, respectively ( j)4.3t . Solutions consistent 
with neutrino radiation transport correspond to the intersection of each curve 
with the T v = 2/3 dashed line. The Burrows-Goshy limit is marked by a thick 
black curve that has a maximum at exactly 2/3. T he ef fective optical depth 
for electron-type neutrinos is used to compute t v ( j|2.3> . 

dynamic point of view^. By relaxing this constraint, one can 
examine the behavior of t v as a function of shock radius and 
neutrino luminosity, in order to gain insight into the processes 
that determine this limit. 

Figure |2] shows the neutrino optical depth r„ as a function 
of shock radius R s at t = for several neutrino luminosities in 
our fiducial sequence. Each point along a given curve is ob- 
tained by integrating the steady-state version of equations (01- 
© from an arbitrary shock position down to a constant R v , 
keeping the upstream flow fixed, and calculating equation ([8]). 
For small shock radii, the optical depth is a monotonically 
increasing function of R s . For large enough shock radii, how- 
ever, t„(R s ) reaches a maximum and then decreases for in- 
creasing R s . Larger neutrino luminosities move this maximum 
towards smaller shock radii and decrease its magnitude. The 
two points where this curve reaches a val ue of 2/3 correspond 
to the two physical solutions found by | Yamasaki & Yamada 
(2005) and Peic ha & Thompson! (1201 ll) . Further increases in 
the neutrino luminosity cause the maximum of r v to fall be- 
low 2/3, making it impossible to satisfy the optical depth con- 
straint. 

The function dr v /dR s has two contributions. The first 
comes from the increase in the volume of the postshock region 
and is always positive for fixed R v . The second arises from a 
change in the quantities that determine k u , chiefly the density 
profile (eqns. I9l- lfl0l ). and can have either sign. Figure [3^ 
shows the density profile of the model at the Burrows-Goshy 
limit and another with neutrino luminosity 5% larger and the 
same shock radius (corresponding to the lowest solid curve in 
Figure[2j. The density of the higher luminosity configuration 
is lower by a factor of the order of 2 near the neutrinosphere, 
where the bulk of the optical depth is generated. 

The density decrease in the cooling layer is caused by an 

4 The settling of a shocked accretion flow onto a spherical surface with an 
arbitrary heating and cooling function is a well-defined hydrodynamic prob- 
lem, albeit not necessarily realistic. 
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FIG . 3 . — Profiles of density (a), ratio of pressure in leptons and radiation to 
total pressure (b), and entropy (c) for two models from our fiducial sequence 
which have the same shock radius but different neutrino luminosities. The 
thick line corresponds to the Burrows-Goshy limit, while the thin line has a 
neutrino luminosity 5% higher. A higher entropy in the gain region increases 
the fraction of the total pressure provided by relativistic particles at fixed 
radius, softening the density profile and thus decreasing the density around 
the neutrinosphere. This is the reason why the lowest curve in Figure EJfalls 
below r„ = 2/3. 

inward motion of the radius at which the pressure transitions 
from being dominated by leptons and radiation to being set 
by non-relativistic nucleons. This is shown by Figure[3]3. The 
adiabatic index in the relativistic particle dominated region 
is closer to 4/3, while nucleon domination moves it towards 
5/3, causing a steepening of the density profile at this tran- 
sitiorfl Inward motion of the transition region relative to R v 
causes the density to be lower at this radius. The change in 
the transition radius in turn is caused by a higher entropy in 
the gain region (Figure [3]:), which increases the contribution 
from relativistic particles. Finally, the higher entropy is the re- 
sult of the higher heating rate arising from the higher neutrino 
luminosity. 

In the case where the neutrinospheric radius is tied 
to the neutrino luminosity at constant temperatu re (e.g., 
iBurrows & Goshylll993L lYamasaki & Yamadall2007l) . an ad- 
ditional contribution arises from the change of the lower limit 
of integration in equation ([8]). This decreases the size of the 
integration volume for increasing neutrino luminosity relative 

5 The steepening of the density profile with a higher adiabatic index occurs 
because the temperature varies slowly with radius near the neutrinosphere. 
The pressure gradient needed to balance gravity thus arises mostly from the 
density gradient. In contrast, in the 7 ~ 4/3 region the pressure gradient 
arises mostly from the changes in the temperature, with the density profile 
being thus shallower I Janka 2001). 



to the case where R v is constant, making t v smaller than it 
w ould otherwise be. 

Peic ha & Thompson! d201 ll) found that the Burrows-Goshy 
limit can be characterized by a critical value of the ratio of 
the sound speed to free-fall speed at some point in the flow, 
above which it is not possible to connect a subsonic settling 
solution to a supersonic accretio n flow via the hydrodynami c 
shock jump conditions (see also Yamasaki & Yamada 2005). 
Implicit in this result, however, is the simultaneous fulfill- 
ment of the boundary conditions at the upstream side of the 
shock and at the neutrinosphere when constructing the two 
solutions. Here we have adopted a uni-directional approach, 
keeping only the upstream boundary condition and integrat- 
ing the fluid equations inward until an additional constraint is 
met at the base. This allows us to violate the optical depth 
condition by varying the shock position, which can be placed 
at any point in the upstream solution. The conclusions of both 
approaches regarding the nature of the Burrows-Goshy limit 
are equivalent when viewed in this light, as a higher sound 
spee d relative to the free-fall velocity implies higher entropy 
(see lMurphv & M eakin 2 01 ll for a discussion of the connec- 
tion between these two quantities through the entropy equa- 
tion). 

4. TRANSITION TO RUNAWAY EXPANSION 

In this section we analyze the instabilities that mediate the 
transition to runaway expansion, using simulations results and 
adapting analysis techniques from stellar pulsation theory. 
The term explosion implies the existence of a star whose en- 
velope is to be ejected by a successful shock. Given the para- 
metric character of our study, we employ instead the term run- 
away expansion to denote the onset of explosion in a medium 
with constant mass accretion rate, as the asymptotic explosion 
energy is determined by addit ional processes not included 
here (e.g. jMarek & Jankall2009h . 

Evolving any of the model sequences described in 32.51 
yields a characteristic outcome for the evolution of the shock 
radius as a function of time. The result for selected models 
from our fiducial sequence is shown in Figure|4] Models with 
low neutrino luminosity are stable to small perturbations. For 
4.05 < L„ f /(10 52 erg s" 1 ) < 4.55, the shock undergoes oscil- 
lations of growing amplitude that transition to runaway ex- 
pansion after some time delay. For L Ve > 4.55 x 10 52 erg s" 1 , 
initial exponential growth in amplitude occurs without any os- 
cillation, transitioning into a nearly constant-velocity expan- 
sion phase at later times. 

The characteristic timescale associated with both types of 
modes is the advection time from the shock to the protoneu- 
tron star 

_ [ R ° dr _M em 

fadv- / 1 1- - — -T-, (ID) 

JR U \Vr\ M 

where the second equality is valid when M is constant, and 
M env is the mass enclosed between the neutrinospheric radius 
and the shock. The fluid velocity is very subsonic, hence the 
kinetic energy content is small. Changes in the accretion flow 
involve changes in the heat content of the fluid through mod- 
ulation of the neutrino emission and absorption processes, 
hence the instabilities involved are non-adiabatic. 
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FIG. 4. — Shock radius as a function of time for selected models from our 
fiducial sequence. Different curves correspond to different electron neutrino 
luminosities as labeled (L Ve j2 =£k,/10 52 erg s -1 ). Stable, oscillatory, and 
non-oscillatory modes are plotted as black, blue, and red curves, respectively. 
The maximum luminosity for a steady-state solution at fixed optical depth 
(the Burrows-Goshy limit, ^3) corresponds to L„ 52 — 4.74 for this set of 
parameters. 

4. 1 . Work Integral for Accretion Shocks 

Any instability involves the conversion of energy from one 
type to another as the system transitions into a more favor- 
able configuration. In stellar pulsation theory, the standard 
tool used to quantify this convers ion is the work integral 
(Eddington ll926t lUnno et afl fl989). which is defined as the 
change in the total energy E of the system over and oscilla- 
tion cycle 

/dE 
At—. (16) 

In pulsating stars, the total mass of the system is essen- 
tially constant within an oscillation period, thus equation ( fT6] l 
is normally wri tten to consider Lagrangian changes (e.g., 
lUnno et alJ[l98M 

In contrast, for an accretion shock neither the mass nor the 
volume enclosed within the star and the shock remain constant 
in time, hence an Eulerian expression must be used if one 
wants to adapt this tool to the problem at hand. Writing the 
total specific energy as 



<?tot = e- 



GM 
r 



(17) 



and rearranging the gravitational contribution to equation ©, 
the Eulerian rate of change of the total energy in spherical 
symmetry becomes 



dE 
~dt 



6?x^^-+AT,R%{pe M )\ R =E tot (18) 



dt 



= E^-E up + Ed n + E s 



with 



En = 



4v?dr% m 



£u P = AttR s [v r (pe tot +p)]\ R 



(19) 

(20) 
(21) 



6 The work integral for stars is also simplified by the pressure falling to 
very small values at the surface, which is not the case for flow confined by an 
accretion shock. 



E dn = 4nR 2 n [v r (pe tot + p)] \ R ^ 
E s =4irR 2 s R s (pe m )\ R , 



(22) 
(23) 



where R s denotes the rate of change of the shock position, and 
Ri n is some inner boundary that remains fixed. Equations ( fT8l - 
( 1231 can be straightforwardly extended to allow for a moving 
inner boundary and non-trivial integrat ion over a ngles. An 
analogous expression was derived by Janka (2001) to study 
the evolution of the mass and energy in the gain region. 

The first three terms in equation (fT9l represent, respec- 
tively, the integrated neutrino source terms plus the energy 
fluxes entering and leaving the integration domain with the 
accretion flow. In steady-state, these three terms cancel out 
exactly. The fourth term is non-zero only when the shock 
moves, and accounts for the change in energy caused by a 
change in the postshock-volume. Because the total specific 
energy just below the shock is generally negative due to dis- 
sociation losses, this term provides damping on expansion and 
driving on contraction. 

Because the evolution of the shocked flow is not necessar- 
ily periodic, we will make use of equation ([Tol l rather than 
( flol i to track the different energetic components in the sys- 
tem. Also, given that the entire flow between R u and R s is in 
sonic contact, analysis of global instabilities must include the 
cooling region, thus we will take as close as possible to 
R v . Details about the calculation of equations (fT8l-(l23l from 
our simulations are provided in Appendix ICl 

4.2. Linear Phase: Oscillatory and Non-Oscillatory Modes 

Figure shows the evolution of the shock radius and the 
different terms that account for change of the total energy 
(eq. |[l9l ). for a model that explodes via an oscillatory mode 
(Lu e ,52 = 4.1). The lower integration boundary R m was cho- 
sen to be a few cells above R v , to eliminate boundary effects 
(Appendix[C]>- 

In the linear phase, the magnitude of the energy genera- 
tion by shock motion E s peaks when the shock velocity is 
the largest, and is dominated by the gravitational term in <? tot . 
Contractions thus inject positive energy into the postshock re- 
gion. The net energy generation from accretion En - E up + £dn 
also becomes non-zero once the system deviates from equilib- 
rium. Its magnitude is much smaller than the individual terms 
that comprise it, showing that imperfect cancellation between 
the energy fluxes and neutrino source terms is the additional 
source of energy. The lowest panel in Figure [5] shows that 
in fact the flow adjusts itself in the linear phase to match the 
energy generation from shock motion with the excess energy 
flux entering the domain, E s (t) ~ E up (t) — E up (0). The total 
energy generation can then be accounted for entirely by the 
deviation from steady-state of the neutrino source terms and 
the energy flux leaving the domain toward the protoneutron 
star. In other words, oscillatory instability arises from a mis- 
match in the dissipation of the accretion energy in a system 
with a moving boundary. 

The origin of this imbalance can be traced back to the 
steep temperature dependence of the neutrino cooling func- 
tion, J^c/p oc T 6 (eq. 1A13I ). Figure|6]shows that the pertur- 
bation to the specific neutrino source term closely mirrors the 
temperature perturbation, with the opposite sign (after remov- 
ing the density dependence, the heating term depends only on 
the neutron and proton abundance, linearly). Thus an increase 
in the temperature from its steady-state value causes an in- 
crease in the rate of cooling and hence excess energy dissipa- 
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FIG. 5. — Evolution of the different rates of energy change (eq. 1191 ) for 
a model in our fiducial sequence which explodes via oscillatory instability 
(L Ue ,52 = 4.1). Panel (a) shows the shock radius as a function of time (solid 
black), and the radius at which the bi ndin g energy of alpha particles equals 
their gravitational binding energy (eq. 1111 . horizontal red dashed). The verti- 
cal dashed black line signals the time at which the fluid achieves positive total 
energy for the first time (see also Fig. |5J, while the green dashed line shows 
the time at which the sound crossing time from the shock to R„ becomes 
longer than the shock expansion time R s / Vj . Panel (b) shows the energy gen- 
eration due to shock motion (blue), the net energy generation from accretion 
(red), and the total rate of change of energy (green). Panel (c) shows the fluc- 
tuation (relative to the initial condition) in the energy flux entering through 
the shock (solid red), the energy generation from shock motion (dashed blue), 
the fluctuation in the energy flux leaving the domain minus neutrino cooling 
(solid green), and the total rate of change of energy (dashed orange). Note 
that the vertical scale changes relative to panel (b). 

tion. This results in contraction and cooling of the postshock 
volume on a sound crossing time, reversing the cycle when 
the lower temperature fluid reaches the cooling region by ad- 
vection. 

Figure [6] also shows that the period of the oscillation is re- 
lated to the advection time from the shock to a radius near the 
point of maximum cooling. At this radius the flow achieves 
maximum deceleration, and advected p erturbations are ef- 
ficiently converted into acoustic waves dSchecket all [2008: 
iFoglizzol 120091; ISato et al.1 12009b . The normalized density 
perturbation is significant only outside of this radius, be- 
low which it undergoes a phase shift. The model shown 
in Figure [6] has a period very close to twice the advec- 
tion timescale from the shock to this radius. This is very 
close to the period of an oscill tory mode with no heating 
dFernandez & Thompson 2009b). As the flow becomes un- 
stable, however, the oscillation period becomes increasingly 
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FIG. 6. — Spacetime diagram showing different perturbed quantities as a 
function of time for a stable model in our fiducial sequence (L u 52 = 3.5). 
Two complete oscillation cycles are shown, and the differences are taken rel- 
ative to values at very late times, where the flow reaches perfect steady-state. 
The top panel shows normalized density fluctuations, the middle normalized 
temperature fluctuations, and the bottom the fluctuation in the specific net 
neutrino source term _S? nc t/p in units of 10 18 erg g s , Horizontal lines cor- 
respond to the gain radius (upper solid white), maximum cooling (dashed 
black/white), and radius where cooling has decreased by an e-folding from 
its maximum (white/black solid near bottom of each panel). Diagonal lines 
correspond to streak lines in the unperturbed solution. The amplitude of the 
shock radius (white wavy curves) has been increased by a factor of 100 for 
clarity. 

longer, achieving nearly four advection times at the threshold 
for oscillatory instability. This lengthening of the oscillation 
period with neut rino luminosity has been a ttributed to buoy- 
ancy effects by Yamasaki & Yamadaj (120071) . 

Non-oscillatory modes grow exponentially in amplitude 
and transition directly into runaway expansion. The differ- 
ent terms that make up the rate of change of total energy are 
shown in Figure|7](the model has L„ f ,52 = 4.6). The evolution 
of these terms resembles the last expansion phase of oscilla- 
tory modes, where the shock motion term acts only as a sink 
of energy and is not matched by the deviation in E U p from 
steady-state. The expansion is driven by the net energy gen- 
eration from accretion, which is positive throughout. 

4.3. Approximate Instability Criteria 

The only way to obtain an exact instability threshold for 
these non-adiabatic modes is to perform a linear pertur- 
bation analysis of e q uation s d4j-(j7j such as that done by 
lYamasaki & Yamadal (120071) . In the absence of such a cal- 
culation, we have instead searched for physically motivated 
instability criteria that correctly describe the behavior of our 
simulations over a wide range of parameter space. We discuss 
here their definition, and leave for ^5] their comparison with 
simulations over an extended region of parameter space. 
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FIG. 7. — Same as Figure \5\ but for a model that explodes via the non- 
oscillatory instability (L„ f 52 = 4.6). Only the two upper panels are shown, as 
the fluctuation in £ up does not compensate E s , as in the last expansion of the 
oscillatory mode. 

For oscillatory modes, we have found that instability sets 
in when the advection time through the gain region is longer 
than the time required to advect from the gain radius to a point 
close to the node in the density perturbation (Figure |6]l times 
some coefficient close to unity. This is equivalent to requir- 
ing more mass to reside in the gain region than in the part 
of the cooling region above the node in the density perturba- 
tion, at any given instant (eq. |[T5l ). To obtain a coefficient of 
unity, the relevant location is the radius at which the density 
and pressure perturbations are 180 degrees out of phase with 
each other. In the steady-state solution, and for most of our 
models, this position is approximately the radius closest to R v 
where the specific cooling has decreased in magnitude by one 
e-folding from its maximum (e.g., FigureQJ). 

The approximate instability criterion is then 



^adv— g ^adv-ei 



with 



and 



^adv— g — 



dr 

H 



^adv-e — / 11 



* dr 

, v > , 

1 ( J2?net 



P 



(24) 
(25) 

(26) 
(27) 



and where the solution to equation (l27l > closest to R u is taken. 
The dependence of these two timescales on neutrino luminos- 
ity is shown in Figure [8] 

It is worth emphasizing that r e was chosen because it yields 
a coefficient of unity. Another location close to the node in 
the density perturbation could also have been used, with a 
slightly different coefficient. For example, taking the radius 
of maximum specific cooling r cmax yields an instability condi- 
tion fadv.gain > 1 -5fcmax, where the latter timescale is computed 




ve, 52 



FIG . 8 . — Characteristic timescales that determine the approximate instabil- 
ity condition for oscillatory and non-oscillatory modes, as a function of elec- 
tron neutrino luminosity in our fiducial sequence (compare wi th F igure fj). 
Shown are the advection time over the gain region (black, eq. 1251 ). advec- 
tion time from the gain radius to the p oint where cooling has decreased by an 
e-folding from its maximum (red, eq. 1261 ), and timescale to c hange the total 
energy in the gain region by neutrino heating (green, eq. 1291 ). The vertical 
dashed line marks the Burrows-Goshy limit. 

by replacing r e with r cmax in equation d26| i. 

A physical justification for r e can be found by looking at 
the behavior of the relative phase between the pressure and 
density perturbations. It is a general result from stellar pulsa- 
tion theory that, for oscillatory instability, an excitation mech- 
anism is required which causes the system to be absorbing 
heat at the point of maximum compression (e.g.. lCoxlll974h . 
This way, the pressure maximum is reached during the expan- 
sion phase, leading to positive work over an oscillation cycle, 
and hence to a net increase in the pulsation kinetic energy. 
In the system under study, such a phase lag of the pressure 
perturbation relative to the density is satisfied below the node 
of the density perturbation. The point where both perturba- 
tions are 180 degrees apart is approximately r e for most of 
our models, with the pressure leading the density at any given 
radius above ifl In cases where r e does not trace the phase 
radius correctly, the criterion in equation d24b loses accuracy. 
The relation between the phase radius and the cooling pro- 
file is most likely a consequence of the way we suppress the 
neutrino source terms with density (eq. lfl4l ): different imple- 
mentations of neutrino transport will most likely find slightly 
different values. 

An interesting question is how the phase lag criterion for 
stars can be adapted to the problem at hand, given that fluid el- 
ements are advected downstream instead of returning to their 
original position. Taken at face value, the phase lag criterion 
would dictate that the driving region resides below the node 
of the density perturbation. However, the pressure lags the 
density in this region for both stable and unstable models, and 
r e hardly changes over a wide range in L Vc (Figure [8]). Such 
question is likely to be at the root of the instability mechanism 
of long waveleng th modes of the Sta nding Accretion Shock 
Instability (SASI. lBlondin et aljr2003l) with heating included, 
and will not be further pursued in this paper. 

In the case of non-oscillatory modes, the instability thresh- 
old is very close to the point where the advection time through 

7 To determine the precedence of density or pressure maxima, we take two 
peaks that are closest in phase. This is only possible when perturbations have 
a phase shift different than or 1 80 degrees. 
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the gain region is longer than the time required to change the 
total energy through heating, 



with 



^adv— g ^ ^heat— g 



eat-g ■ 



J R ^d 3 x(pe tot ) 



(28) 



(29) 



The dependence of fheat-g on neutrino luminosity is also shown 
in Figure [8] This relation has long been known to provide 
a predictive criterion for runaway e xpansion in spherically 
symmetric core-collapse simulations (iThompson et alj|2005t 



Buras et al. 2006a; Murphy & Burro ws 2008; Mar ek & Jankal 
20091) . Itsinterpretation is straightforward: neutrino heating 



is effective enough to change the internal energy of the flow 
during its transit time through the gain region. As a conse- 
quence, the pressure cha nges significantly an d the shock re- 
adjusts to a new position ( Janka & Keil 1998|). 

Note that equation d28l i is a global condition on the gain 
region and includes the gravitational binding energy, dif- 
fering from local defi nitio ns of this ratio such as tho se in 
Thom pson et al.1 (120051) and lPejcha & Thompson! (1201 ll) . The 
fluid does not achieve positive energy over a time f at jv-g (~ 
10 ms, Figure|8}, but instead takes a multiple of this timescale 
to transition into runaway expansion. As we discuss in the 
next section, onset of runaway occurs when the fluid achieves 
positive energy for the first time. 

We have checked whether any instability threshold can be 
described by a fixed value of the parameter \ that measures 
the effects of buoyancy in an advective flow, 



* dr 

. Wo 



where Wbv is the Brunt- Vaisala frequency, 



2 _ GM ( 1 ainp 



dlnp 

dr 



(30) 



(31) 



with T s = {p/p)c 2 s . When \ > 3, buoyancy is expected to over- 
come the stabilizing effect of advection and l ead to connec- 
tive i nstability in the multidimensional case (Fogliz zo et al.l 
2006). Even though all of our sequences have \ in the range 
1-10 near the threshold for instability, in none of them does 
this parameter have a constant value along the critical stabil- 
ity curves. For example, in our fiducial sequence, the value 
of this parameter at the non-oscillatory threshold ranges from 
X~7 forM = 0.1M Q s" 1 to \ = 1 .7 atM = 1M W s" 1 . It is worth 
pointing out that the analysis of Foglizz o" et al.l (12006b applies 
to infi nitesimal perturbations. The results of IScheck et al.l 
(2008) show that the SASI can trigger convection through fi- 
nite amplitude density fluctuations in cases where x < 3, so 
our findings do not necessarily imply that convection will be 
suppressed at large accretion rates in the quasi steady-state 
approximation (§2. 1). 

4.4. Nonlinear Phase and Runaway Expansion 

Once oscillatory modes become unstable, the amplitude 
grows steadily until the oscillation cycle is broken and the sys- 
tem transitions into runaway expansion. Figure [5] shows that 
not only the shock radius but also the different energy genera- 
tion terms undergo a qualitative change in behavior relative to 
the linear phase once this point is reached. In particular, the 
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FIG. 9. — Total specific energy e to t (a) and radial velocity (b) as a function 
of radius for the model shown in Figure[5] which explodes via the oscillatory 
instability. Different curves correspond to different times in the simulation, as 
labeled. The black curve (277 ms), red curve (292 ms), and green curve (314 
ms) correspond to the vertical dashed lines of the same color in Figure[5ji. 

energy generation due to shock motion decouples from the 
fluctuation in the energy flux entering though the shock. Ex- 
pansion is driven by the net energy generation from accretion, 
with damping from the shock motion term. 

A physical origin for this transition point can be found by 
inspecting the evolution of the total specific energy and radial 
velocity in the system. Figure [9] shows that this time corre- 
sponds approximately to the point where the specific energy 
behind the shock becomes positive for the first time. This 
time is marked by a vertical black dashed line in Figure [5^ 
and the snapshot at 277 ms in Figure [9] At about the same 
time, the fluid just inside the shock achieves positive radial 
velocity, reversing the accretion flow. This transition point is 
equivalent to the esc ape temperature condition discussed in 
iBurrowset all (119951) . 

Physically, the above result is dependent on the definition 
of the zero points of energy in the context of a Newtonian 
framework. The EOS used her e takes this zero p oint to be, per 
baryon, the atomic mass unit (She n et al.lll998l) . Our expres- 
sion for the total specific energy, equation ( fT71 ), takes the zero 
point of gravitational binding energy at an infinite distance 
from the central mass. Given these definitions, the condition 
of positive energy leading to runaway expansion is a well- 
defined mathematical concept. A self-consistent determina- 
tion of this condition would require inclusion of the rest mass 
energy and gravitational field calculated self-consistently in a 
general relativistic framework. 

The shock accelerates when it approaches the radius where 
the binding energy of alpha particles equals their gravitational 
binding energy (eq. ifTTI ). This can be understood from the 
fact that the dissipation due to shock motion E s decreases in 
magnitude with increasing radius, because the total specific 
energy behind the shock becomes less negative. This can be 
seen in Figure [5}), where the blue curve reaches a minimum 
value. By now most of the gain region has reached positive 
energy, except for a narrow layer behind the shock, and most 
of the fluid has positive velocity, effectively becoming a wind- 
like solution. 

A final transition occurs when the shock begins to mechan- 
ically decouple from the protoneutron star atmosphere. This 
occurs approximately at a time when the sound crossing time 
from the shock to R v becomes longer than the shock expan- 
sion time R s /v s , where v s is the shock velocity. This time is 
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FIG. 10. — Same as Figure [9] but for the model shown in Figure [7] which 
explodes via the non-oscillatory instability. Black (64.5 ms), red (75 ms), 
and green (92 ms) curves correspond to the times signaled by vertical dashed 
curves of the same color in Figure|7] 

marked as a vertical green dashed line in Figure [5^ and cor- 
responds to the curve at 314 ms in Figure|9] This mechanical 
decoupling explains the fact that the net energy generation 
becomes negative yet the shock continues to expand, as in- 
ferred from Figure[5] It is worth noting that velocities below 
the shock are subsonic throughout the time period considered 
here. 

Non-oscillatory modes transition directly into runaway ex- 
pansion, and their evolution closely resembles that of oscilla- 
tory modes after reaching positive energy. The corresponding 
snapshots of total specific energy and radial velocity at differ- 
ent phases are shown in Figure flOl 

In all of the sequences simulated, transition to runaway ex- 
pansion is achieved after the flow becomes radially unstable. 
We have found no unstable mode that saturates, nor any model 
that having achieved positive energy, failed to continue into 
runaway expansion. Our preliminary conclusion is that for 
the set of assumptions adopted in this paper, radial instability 
is a sufficient condition for transitioning into runaway expan- 
sion. 

4.5. Conditions for Saturation 

There are many examples in the literature where one- 
dimensional stalled shocks start to expand or os cillate, but 
then fizzle (e.g., Figure 1 of|Janka & Miiller 1996). It is thus 
imperative to identify the processes that can lead to saturation, 
as our preliminary finding linking instability to runaway may 
break down in some region of parameter space. 

Insight on this question can be obtained again by inspecting 
the different processes that contribute to the change in the total 
energy in the post-shock region (eq. fT9l ). What we are inter- 
ested in are dissipation processes that contribute with negative 
energy generation. 

In the present study, we are ignoring evolution of the neutri- 
nospheric parameters (R v , core luminosities, and spectra) and 
heating due to the accretion luminosity, hence the dominant 
dissipation mechanism arises from the change in the post- 
shock volume, E s , as inferred from Figures[5]and[7] This term 
is made up of three factors: (1) the specific energy below the 
shock, (2) the density profile, and (3) the rate of change of the 
post-shock volume. 

1. At the typical radii where supernova shocks stall (100- 
200 km), the total specific energy <? tot is usually nega- 



tive. Because the temperature decreases outward, nu- 
cleons below the shock recombine first into alpha par- 
ticles and then heavy nuclei as the shock moves out, 
yielding an increase in the thermal energy t hat is com- 
parable to the gravitational binding energy (BetheJUH 
Fern andez & Th ompson 2009a). This causes the total 
specific energy to become significantly less negative, 
decreasing the size of E s as the shock expands. This can 
be most clearly seen in Figure \5&, which shows a clear 
transition in the evol ution of E x when the shock reaches 
r a . The results of iFernandez & Thompson! ((2009a) 
show that using a constant dissociation energy indeed 
quenches runaway expansion, in direct contrast to al- 
lowing alpha particles to recombine, because in the first 
case the dissociation energy becomes an increasingly 
larger fraction of the local gravitational binding energy 
as the shock expands. Even then, saturation occurs at 
a radius that is several times the initial shock radius. 
Except for unrealistically small shock stagnation radii 
(< 50 km), this saturation channel is unlikely to be of 
importance. 

2. The evolution of the density profile below the shock 
and the mass accretion rate are tied to the density pro- 
file of the progenitor at the onset of collapse. For an 
iron core supported by relativistic electrons the den- 
sity scales like r~ 3 , yielding a mass accretion rate that 
scales inversely with time at fixed radius when mass 
conserv ation and a n ear free-fall velocity field are as- 
sumed dBethd [T990I) . A time-independent mass ac- 
cretion rate corresponds to a progenitor density profile 
oc r~ 3 / 2 , which would be given by an adiabatic index 
of 5/3 in hydrostatic equilibrium. In other words, a 
time-independent mass accretion rate, as assumed in 
our sequences, is already unrealistically high and over- 
estimates the dissipation. 

3. Because E„ is negative on expansion, it tends to sta- 
bilize the shock velocity for a given energy dissipa- 
tion rate. The dominant energy source for expansion is 
neutrino heating. As long as the total heating remains 
larger than E s before the shock is completely decoupled 
from the atmosphere, the shock will continue to expand. 

We therefore conclude that for the physical assumptions of 
the present study, radial instability is a sufficient condition for 
runaway expansion. 

In the general case, however, evolution of the neutri- 
nospheric parameters and neutrino heating by the accretion 
luminosity can provide significant energy dissipation and sat- 
urate the instability. By inspection of equation (TE[ , one can 
find three sources of negative energy. 

First, evolution of the core neutrino flux affects the energet- 
ics via £n- If the net effect of decreasing neutrino luminosities 
and increasing mean neutrino energies is a decrease in the en- 
ergy deposition in the gain region, the instability can be sup- 
pressed in two ways. The linear instability depends on the ex- 
tent of the gain region through ? a dv,gain, hence a rapid decrease 
of heating can stabilize oscillations or non-oscillatory expan- 
sion if the flow has not yet achieved positive energy when the 
stability criterion is reversed. Second, a decrease in the rate of 
neutrino energy deposition can lead to quenching of the run- 
away phase if the total rate of heating fails to keep up with 
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FIG. 11. — Threshold electron neutrino luminosity for radial instability as a function of mass accretion rate, for different sequences of models. Panels (a), 
(b), and (c) correspond to sequences with fixed neutr inospheric radii at R v = {30,20,40} km, respectively, while the sequence in panel (d) relates R v to the 
neutrino luminosity via the black body relation ( 82.51 . Data points denote results from time-dependent simulations, where the threshold for explosion is taken 
as the average between stable and unstable models, with the error bars showing the size of the difference in L Vr (AL„ e < 10 51 erg s ). Blue points c orre spond 
to oscillatory and red to non-os cilla tory instability, respectively. Curves correspond to our approximate instability criteria for oscillatory (blue, eq. 1241 ), and 
non-oscillatory modes (red, eq. 1281 ). Black curves correspond to the limiting luminosity of the steady-state configuration at constant optical depth, calculated as 
described in ij3] Black circles in panel (b) denote the absence of instability (and thus explosion) for M > 0.5Mg s for all luminosities up to and including the 
Burrows-Goshy limit. 

dissipation due to E s . The results of Uanka & Mullerl (119961) 
are consistent with the operation of this saturation channel. 

Second, the contribution of the accretion luminosity to 
the total neutrino flux i n realistic models is ~ 50% (e.g., 
Liebendorferet al. 2001). The important aspect to keep in 
mind here is that shock expansion relative to its steady-state 
position leads to a decrease in the magnitude of the mass ac- 
cretion rate. Perturbin g the mass conservation equation yields 
dFoglizzo et alj|2007l) 



try. A larger core neutrino flux is therefore required to sus- 
tain expansion, relative to the light bulb heating case, and the 
causality relation between radial instability and explosion we 
have found here will be violated. 

Finally, the contraction of the protoneutron star generates 
in itself a sink of energy in the post-shock region. Including 
this effect, and keeping everything else constant, would add a 
term of the form 



SM 
M 



1 



8v, 



(32) 



-E^sph = -AirR u R v (pe lot ) R 

"\50n2 



(33) 

1 0"7?f y M R vfi puew.w ergs" 1 (34) 



where vi and V2 are the (negative) upstream and down- 
stream fluid velocities, and Sv s is the shock velocity pertur- 
bation. Taking Sv s real and positive (shock expansion) yields 
SM/M < 0, or a decrease in the magnitude of the accre- 
tion rate. Hence including the heating from accretion neu- 
trinos adds a non-trivial feedback to both oscillato ry and non- 
oscillat ory modes, altering the instability criteria. Bur as et alj 
(2006b) find that indeed spherical oscillations are damped by 
the decrease in the heating due to the dropping mass accretion 
rate. As noted by a number of previous works, the runaway 
expansion phase also cuts off accretion in spherical symme- 



to equation ( [T8l ). with R m = R„. The notation in the sec- 
ond equality is R v3(l = R u /(30 km), R vA = R v /(-\Q b cm s _1 ), 
pu = p/(10 n g cm -3 ), and e to t,i9 = e to t/(-10 19 erg g _1 ). In 
contrast to E s , contraction leads to energy dissipation. This 
quantity is smaller than typical core neutrino luminosities, 
but approache s the net energy generation terms that regulate 
instabilities ( 34.21 ). A more careful analysis would need to 
include the additional cooling and heating by accretion neu- 
trinos generated by enl arging the post-shock c avity. The one- 
dimensional results of Janka & Miiller ( 1996|) show that in- 
deed ~ 10% higher core neutrino luminosities are required to 
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start an explosion in models that experience protoneutron star 
contraction relative to fixed cores, when all other parameters 
are similar. 

5. INSTABILITY THRESHOLDS AND RELATION TO THE LIMITING 
STEADY-STATE LUMINOSITY 

5.1. Dependence on Mass Accretion Rate and 
Neutrinospheric Radius 

For each of the simulation sequences described in 32.51 we 
have searched for the instability thresholds of oscillatory and 
non-oscillatory modes. Figure [TT| shows the resulting thresh- 
old luminosities as a function of mass accretion rate. In all 
cases, data points correspond to the average luminosity be- 
tween two models at both sides of the instability threshold. 
The separation in luminosity, shown as error bars, is less than 
10 51 ergs _I . 

Figure QT] also shows the limiting luminosity for a steady- 
state configuration (the Burrows-Goshy limit), calculated as 
described in $3] The resulting curve agrees qualitatively with 
the results oflBurrows & Goshvld 19931). [Yamasaki & Yamada 
(120051 12006b . and lPeicha & Thompson] d201 lb . Note that the 
neutrinospheric temperatures, neutrino opacities, and EOS 
employed here differ from what was used in those studies, 
hence numerical values are expected to differ. 

The calculation of this limiting luminosity in a consistent 
manner with the microphysics, initial conditions, and bound- 
ary conditions employed in our simulations allows direct test- 
ing of the hypothesis that the critical stability threshold for 
explosion is given by the Burrows-Goshy limit. We find that 
in all cases, the measured thresholds for both type of mode lie 
below this limiting value. When the neutrinospheric radius is 
held constant, the instability thresholds approach (but never 
equal) the Burrows-Goshy limit for increasing mass accretion 
rate. In the sequence that relates the neutrinospheric radius 
to the neutrino luminosity via the black body relation (Fig- 
urefTTtl). the instability thresholds have a nearly constant sep- 
aration in luminosity over the entire range of accretion rates 
investigated. 

The sequence with R u = 20 km yielded unexpected results 
for M > 0.5M© s" 1 , however. In this case there is neither in- 
stability nor runaway for all luminosities up to and including 
the limiting value, independent of numerical resolution. Our 
approximate instability criteria fail to predict this behavior, 
showing that additional constraints play a role in determin- 
ing instability. Also, the fact that no explosion is found at 
the Burrows-Goshy limit shows that this luminosity is not an 
independent tracer of stability either. Increasing the neutrino 
luminosity above the limiting value (at constant shock radius) 
by ~ 10% causes the shock to readjust to a new equilibrium 
position, which is also stable. It is worth noting however that 
the shock radius is R s < 61 km for the non-exploding seg- 
ment of this sequence. This is an unrealistically low value, 
with advection times of the order of ~ 3 ms, a factor of a few 
from the dynamical time. It is like ly th at the fixed value of r v f 
used to set the upstream velocity ( 32.31 ) leads to these extreme 
conditions. 

In all other sequences, the approximate instability criteria 
found in 34.3| provide a good description of the measured val- 
ues over a wide region of parameter space, with agreement 
better than 5% in L Ve . They correctly capture the disappear- 
ance of the oscillatory mode in the R u = 30 km sequence at 
high accretion rates, where only a non-oscillatory mode is 
measured. This phenomenon occurs because at large accre- 



tion rates, the shock radius becomes increasingly smaller. The 
size of the resulting gain region is such that the advection time 
through it never exceeds f a dv-e before the heating time drops 
below f a dv-g due to the increasing luminosity. The deviation of 
the measured oscillatory threshold from the approximate cri- 
terion at high accretion rate in the R v = 40 km sequence is due 
to r e not being a good tracer of the point where pres sure and 
density perturbations are 180 degrees o ut-of-phas e ( 34.31 ). 

Our results extend the findings of Yamasaki & Yamada 
( 120071) . who obtained £ = instability thresholds below the 
Burrows-Goshy limit, to a wider region of parameter space. 
Our results are quantitatively different from theirs, partly 
due to the different neutrinospheric temperatures that they 
employed (T Ve = Tp e = 4.5 MeV), their use of the pure ab- 
sorption c oeffi cient for computation of the optical depth (see 
32.31 and 35.21 ), and their inclusion of self-gravity . Using the 
same parameters as Yamasaki & Yamada (2007), we obtain 
a limiting luminosity of £bg,i/ c = 7.36 x 10 52 erg s" 1 for M = 
IMq s" 1 , which differs by a factor of nearly two, and is close r 
to what was found originally by iBurrows & Goshvl dl993l) . 
The corresponding approximate instability criteria for oscilla- 
tory and non-oscillatory modes are 6.41 x x 10 52 erg s" 1 and 
7.04x xl0 52 ergs _1 . These values are lower than our limiting 
luminosity by 13% and 5%, respectively. 

5.2. Other Parameter Dependencies 

Given the parametric character of this study, certain choices 
had to be made in order to construct a background flow that 
is as realistic as possible. Here we explore how our results 
depend on numerical resolution, boundary conditions, param- 
eters that determine the upstream flow, and the suppression of 
source terms near the neutrinosphere. 

Figure [12k shows the approximate instability criteria and 
Burrows-Goshy limit for our fiducial sequence, together with 
instability thresholds measured from simulations at different 
resolutions. The grid is chosen logarithmically spaced, with 
a ratio of spacing between adjacent cells £ = {Rmzx/Rv) 1 ^', 
and Ar m ; n = R V {C,~ 1)> w i m Ar m ; n the cell adjacent to the in- 
ner boundary (Appendix |B]). Oscillatory modes asymptote to 
the theoretical threshold for increasing resolution, albeit con- 
vergence is non-mono tonic. Non-oscillatory modes converge 
monotonically to a slightly different value, indicating that the 
instability criterion is indeed approximate. 

The results shown in Figure Q~2] also indicate that pre- 
vious hydrodynamic studies that used similar methods to 
solve the hydrodynamic equations {Murphy & Burrowsl2008l : 
iNordhaus et al.ll2010UHanke et al.11201 11) have enough resolu- 
tion to capture the critical stability thresholds correctly. How- 
ever, we have found in our models that there are noticeable 
differences in the growth rates of oscillatory modes as a func- 
tion of resolution. For the model with L Vt $x = 4. 1 in our fidu- 
cial sequence, the time delay to explosion varies by ~ 100 ms 
when going from N r = 240 to N r = 3200. True convergence, in 
the sense that the fractional deviations in the shock radius evo- 
lution are much less than unity for models at different resolu- 
tions, is only achieved for N r > 1600 with our fiducial param- 
eters. It is not straightforward to prescribe a definite resolu- 
tion required for convergence given a hydrodynamic method, 
however, as the implementation of the microphysics (tabular 
in our case) also influences this number. A set of standard hy- 
drodynamic tests tailored specifically for the supernova prob- 
lem would be helpful in assessing the reliability of results as- 
sociated with a given implementation. 
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FIG. 12. — Dependence of the approximate instability thresholds and 
Burrows-Goshy limit on various parameters, taking our fiducial sequence as 
baseline. Solid line s denote the approximate instability co ndit ion for oscilla- 
tory (blue, eq. 1241 ) and non-oscillatory modes (red, eq. (28)), respectively, 
while the black solid line denotes the Burrows-Goshy limit. Panel (a) shows 
simulation results (as data points) as a function of radial resolution in a log- 
arithmic grid from 30 km to 1000 km. Other panels show the dependence of 
these critical points on the upstream electron fraction (b), upstream entropy 
(c), radius at which upstream veloc ity i s set to the free-fall speed (d), and den- 
sity cutoff for source terms (d, eq. 1141 ). Parameters of our fiducial sequence 
are denoted by a black dashed line. 

To test the influence of the steady-state boundary condition 
on our results, we have evolved a sequence with a different 
inner boundary condition, which allows an arbitrary amount 
of mass to leave the domain, while providing enough pres- 
s ure s upport to prevent the shocked envelope from collapsing 
( j32.4l >. Figure[T3h shows that the system still goes through the 
same stability phases as the neutrino luminosity is increased, 
with a small quantitative difference < 5% in the critical stabil- 
ity points relative to our fiducial sequence. We thus conclude 
that the two types of instability do not depend on the fluxes of 
mass, momentum, and energy leaving the domain being fixed. 
The quantitative differences can be explained by the fact that 
the outflow boundary condition acts as a persistent source of 
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FIG. 13. — Evolution of the shock radius in selected models of our fidu- 
cial sequence using different prescriptions for the inner boundary condition. 
Panel (a) shows models with dif fere n t ne utrino luminosities using the out- 
flow boundary condition (eqns. 1121 - 1131 ). The curves shown bracket the 
instability threshold for oscillatory (black and solid blue) and non-oscillatory 
(dashed blue and red) modes. Except for the boundary condition, parameters 
are identical to those in Figuref4] Panel (b) compares two runs with the same 
parameters except the inner boundary condition. Both models are stable and 
do not explode (note the vertical scale). 



waves, keeping shock oscillations at some non-zero ampli- 
tude. This effect is shown in Figure[T3b. which compares the 
shock radius evolution for two stable runs with identical pa- 
rameters except for the boundary condition. 

This small quantitative difference, caused by a persistent 
wave source, is interesting as it provides an illustration of 
what we expect will occur in the multidimensional case. The 
critical stability thresholds are modified in the presence of 
this additional wave generation, as comparison of Figures [4] 
and [13] show. Taking the time-average of two stable runs 
with differing boundary conditions after they have settled into 
steady- or quasi-steady-state, shows that the mean velocity 
profiles differ (Figure fl4l>. The magnitude of the change in 
the threshold luminosities is consistent with the magnitude 
and direction of the changes in the approximate instability 
criteria, using the mean advection times rather than the un- 
perturbed ones. A stronger source of wave energy (e.g., con- 
vection) will affect the mean advection time in stronger way, 
and critical stability points will move accordingly. 

Figures [TZb-e illustrate the sensitivity of the approximate 
instability thresholds and t he B urrows-Goshy limit to the pa- 
rameter choices made in 32.31 Shown is the dependence on 
the upstream electron fraction, upstream entropy, radius at 
which the upstream flow has free-fall speed, and cutoff den- 
sity for neutrino source terms (eq. iPBlD . There is a somewhat 
sensitive dependence on the upstream entropy, which in a re- 
alistic situations will be set by the entropy in the progenitor 
core. Note that larger luminosities are needed to cause stalled 
shocks to become unstable when the progenitor entropy is 
smaller. Nonetheless, the relation between the approximate 



HYDRODYNAMICS OF CORE-COLLAPSE SUPERNOVAE 



15 




-10 1 < 1 < < < J — 1 

30 50 100 

radius [km] 

FIG. 14. — Time averaged radial velocity as a function of radius for two 
stable models in our fiducial sequence with the same parameters except the 
inner boundary condition. The time average is taken after the models have 
settled into a steady- or quasi-steady state. 



thresholds and the Burrows-Goshy limit persists without qual- 
itative changes. All other parameters produce changes smaller 
than ~ 5% over a large range of values. 

We have also examined the robustness of the hierarchy be- 
tween instability thresholds and Burrows-Goshy limit when 
the quantities used to calculate the latter are changed. Fig- 
ure [15^ shows the effect of increasing or decreasing the neu- 
trinospheric optical depth (eq. [8)) by a factor of two. Neither 
of these changes cause the Burrows-Goshy limit or instabil- 
ity thresholds to cross each other, for a given optical depth. 
This result can also be seen from the non-overlapping curves 
in Figure [2] for fixed accretion rate. Interestingly, the exact 
optical depth chosen seems to matter the least when going to 
low mass accretion rates. 

The dependence on the choice of neutrino opacity is shown 
in Figure [15b . Here we compare the effects of using the 
pure absorptio n, e ffective, or total (absorption plus scatter- 
ing) opacity ( jj2.3l ) in equation ([8]). Again, the three curves 
maintain their relative hierarchy and do not cross for a given 
choice of opacity. 



5.3. Relation between the Burrows-Goshy Limit and the 
Instability Thresholds 

We have found that the difference between the Burrows- 
Goshy limit and the radial instability thresholds is finite and 
measurable, and is independent of parameter choices. How- 
ever, the question remains as of why is this limiting luminosity 
always close to the threshold for non-oscillatory instability in 
the parameter region relevant to core-collapse supernovae. 

In a more e xtended parameter space st udy of steady-state 
configurations, Pejcha & Thompson] (1201 lh found that at the 
Burrows-Goshy limit, the ratio between f a dv-g and the global 
heating timescale (without gravity) is always close to unity. 
Ultimately, this proximity must depend on the characteristic 
magnitude of the neutrino opacities, the gravitational binding 
energy around a forming neutron star, and the fact that the 
flow is supported against gravity by the thermodynamic pres- 
sure of the gas (e.g., not centrifugal forces). In other words, 
it could be just a result of dimensionality. A possible way to 
test this hypothesis is to explore the behavior of the system 
at lower mass accretion rates, where the expected trend is an 
increasing deviation. 
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FIG. 15. — Approximate instability thresholds and Burrows-Goshy limit 
as a function of mass accretion rate for different optical depths (top) and 
different prescriptions for the neutrino opacity (bottom). The three curves 
remain separated. 

6. SUMMARY AND DISCUSSION 

We have investigated the transition to runaway expansion 
of a stalled core-collapse supernova shock in spherical sym- 
metry, when the parameters that describe the accretion flow 
are varied systematically. A realistic equation of state and 
weak interactions were employed to perform time-dependent 
simulations with FLASH3.2. Starting from steady-state solu- 
tions to the hydrodynamic equations, we evolve sequences of 
time-dependent models with increasing neutrino luminosity, 
and analyze the hydrodynamic processes that mediate the 
transition from accretion to runaway expansion. Our findings 
can be summarized as follows: 

1. - The onset of radial instability is a sufficient condition for 
runaway expansion when heating by the accretion luminosity 
is ignored and neutrinospheric parameters remain constant 
in time. Radial instability can manifest itself via oscillatory 
and non-osci llatory modes, as fou n d in th e linear stability 
analysis of Yamasaki & Yamada (120071) and numerous 
time-dependent hydrodynamic studies. These modes are 
non-adiab atic, as they involve changes in the heat content of 
the fluid ( jg2) . 

2. - For both types of modes, transition to runaway expansion 
occurs after a portion of the fluid in the gain region achieves 
positive energy. This coincides with the fluid just below 
the shock achieving posi tive velocity, starting the phase of 
runaway expansion ( jj4.41 Fi gures l9l and ITOb . 
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3. - The only significant source of dissipation in our models 
is the energy loss from shock motion (eq. |23~1 : Figures and 
[7]), which provides damping on expansion. Nuclear recom- 
bination and a steady source of neutrino heating combine 
to ensure that this term never extinguishes the runaway or 
saturate the linear instability. In a more realistic context, the 
dominant dissipative processes are the decrease of the core 
neutrino luminosity with time, and self consistent heating by 
the accretion luminosity. The latter is expected to provide 
a negative feedback during ex pansion due to the de crease 
of the mass accretion rate (e.g., Janka & Miiller 1996). The 
contraction of the protoneutron star is a somewhat smaller 
correction, of the orde r of ~ 10%. None of these effects is 
included in this study Q4.51 >. 

4. - We have found approximate instability criteria for 
oscillatory and non-oscillatory modes that correctly describe 
the behavior of the system over a wide region of parameter 
s pace , with a precision better than 5% in neutrino luminosity 
( j34.3t Figure ITTb . For oscillatory modes, instability arises 
when the advection time over the gain region becomes longer 
than the advection time from the gain radius to the point 
where the pressure and density perturbations are 180 degrees 
out of phase (Figures [6] and [8). In our implementation, 
this point lies at a position in the background flow where 
the cooling has decreased by an e-folding from its peak 
value. This equality of advection times is equivalent to 
an equality of masses of the respective regions when the 
accretion rate is constant. We have not been able to identify 
a conclusive physical reason for why this condition on the 
masses or advection times triggers an oscillatory instability. 
We surmise that a relation might exist between the fraction 
of the oscillation cycle that the fluid gains ene rgy and th e 
phase lag criterion in stellar pulsation theory (e.g.. lCoxll974l) . 

4. - Non-oscillatory modes become unstable when the ad- 
vection time through the gain region becomes longer than the 
integrated total energy divided by the integrated net heating in 
the gain region. This condition means that heating is effective 
at increasing the thermal energy while the fluid transits 
the gain region, increasing the pressure an d causing the 
shock to adjust to a new equilibrium position dJanka & Keill 
1998). This global criterion has been used b y a number 
of previous studies as an explosion diagnostic (jBuras et alj 
20061 iMurphv & Burrowsll2008t iMarek & Jankafl2009l) . 

5. - The instability thresholds are in general different 
from the limiting lumin osity for the steady-state system 
dBurrows & Goshy||1993L §|3). For constant neutrinospheric 
radius and increasing accretion rates, the thresholds asymp- 
tote to the limiting luminosity, but never coincide with it 
(Figure fTTT) . This separation does not depend on how this 
limiting luminosity is calculated (Figure ITSb . or on specific 
parameter choices (Figure fT2l. 

6. - The existence of a limiting luminosity for steady-state 
solutions is a direct consequence of requiring the optical 
depth between the neutrinosphere and the shock to have a 
fixed value (Figure 13. For fixed upstream conditions and 
increasing neutrino luminosity, the entropy increases and 
hence the radius at which the pressure transitions from being 
dominated by relativistic particles to being dominated by 
non-relativistic nucleons moves inward (Figure 0). This 
causes the density profile to soften at fixed radius, resulting 



in a lower density at the neutrinosphere. Above a certain 
limit, there is not enough mass to provide sufficient neu- 
trino optical depth to satisfy the closure relation, and no 
steady- state solution is poss i ble. T his result is equivalent to 
that of iPejcha & Thompson! (1201 ll) . differing only in which 
boundary conditions are assumed to be fulfilled. 

7. - We find neither instability nor explosion in our sequence 
with constant R u = 20 km forM > 0.5M Q s" 1 and luminosities 
up to and including the limit luminosity (Figure [Tib). At the 
lower end of these mass accretion rates, the shock radius is 
R s ~ 60 km, and decreases for larger accretion rates. We did 
not find an explanation for this result within our framework. 
Regardless of the reason, however, it shows that the transition 
to runaway expansion does not necessarily occur at the 
limiting luminosity. It also shows that the approximate 
instability criteria derived here are incomplete, as an addi- 
tional constraint is likely to determine the stability of the flow. 

Our results show that the lBurrows & Goshvl (11993b conjec- 
ture, relating the transition to explosion to a global instability 
of the shocked envelope, is correct within a restricted set of 
assumptions. The critical stability surface, however, is not 
given by the limiting luminosity of the steady-state configu- 
ration. This limiting luminosity remains nonetheless close to 
the instability threshold over a large section of the parameter 
space relevant to core-collapse supernovae, a result that we 
have not accounted for and which deserves further investiga- 
tion; 

IMurphv & Burrowsl ((2008) and Fernan dez & Thompson! 
(2009a) have also found oscillations around the transition to 
explosion. However, the width in neutrino luminosity separat- 
ing exploding from non-ex ploding configurations is larger i n 
those studies. In the case of Fernandez & Thompson (2009a), 
this can be attributed to the differences in the employed micro- 
physics relative to the present study. La rge amplit ude oscil- 
lations are excited in the simulations of Murphy & Burrowsl 
(2008) when the Si/O composition interface is accreted 
through the shock. Stable models damp oscillations, and un- 
stable modes explode, with neutrally stable oscillations con- 
fined to a range AL„ f jL Vc ~ 4%. This range is wider than that 
found here, presumably because Mur phy & Burrowsl (120081) 
include the contraction of the protoneu tron star. Ou r result s 
are in good agreement with those of Ohn ishi et al.l (12006ft . 
who used essentially the same physical assumptions as we do. 

The main result of this paper, point (1) above, does not ap- 
ply to all modes of higher dimensionality. Non-radial oscil- 
latory instability of the shock (the SASI) does not necessarily 
lead to explosion when the same set of ass umptions adopted 
in this study are employed (Ohnis hi et al.l 120 06). Given our 
results and those of lYamasaki & Ya mada ( 2007), the question 
then arises as to whether a multi-dimensional explosion can be 
thought of as the excitation of an unstable radial-oscillatory, 
and/or non-oscillatory mode (radial or otherwise) by the ac- 
tion of turbulent stresses from the SASI and convection. In 
this case, the excited mode would arise from a background 
state that differs from the laminar steady-state solution by the 
presence of a turbulent pressure term in the momentum equa- 
tion and a convective flux term in the energy equation. The 
excited mode would also become unstable at a lower neutrino 
luminosity than in the spherically symmetric flow. 

According to Yam asaki & Yamadal d2007l) . the non-radial- 
non-oscillatory modes that have the lowest threshold luminos- 



HYDRODYNAMICS OF CORE-COLLAPSE SUPERNOVAE 



17 



ity for instability have Legendre indices I ~ 6. These modes 
are likely to be associated with convection in the gain re- 
gion, as they have no unstabl e osc illatory counterpart. T he re- 
sults of lOhnishi et al.1 (120061) and llwakami et all (12008b show 
that small-scale convection in itself does not lead to runaway 
expansion. In contrast, an I = 1 or £ = 2 non- oscillatory 
mode, such as that envisioned by Thompson (2000), is likely 
to be behind unipolar or bipolar explosions seen in axisym- 
metric simulations (e . g.. IScheck et aTll2006l) . The results of 
lYamasaki & Yamadal ( 120071) show that the threshold lumi- 
nosities of these large scale non-oscillatory modes are larger 
than that of I > 3 modes, and very close to that of the spherical 
oscillatory mode. 

Which of these modes is excited first in a multidimensional 
context will obviously depend on the nature of the new back- 
ground flow. The e xploding three-dimension al models of 
iNordhaus etal] (120101) and IHanke et al.1 (120111) lack a dom- 
inant oscillatory mode in their transition to explosion, even 
though oscillations in the average shock position are still vis- 
ible (particularly for the L Ve = 9 x 10 51 erg s" 1 m odel in the 
11.2M Q progenitor sequence of IHanke et al.ll20TTb . In con- 
trast, two-dimensional models have a noticeable radial oscil- 
latory component, with successive dips in the average shock 
radius before and during runaway expansion. Yet the ampli- 
tude of these oscillations is smaller than in the spherically 
symmetric case, suggesting that more than one mode is likely 
to be involved. 

Support for this interpretation of multidimensional ex plo- 
sions can be found in the results of B uras et alj d2006al) and 
iMarek &~a nka (2009), who use the same definition of fh ea t-g 
as we do here. They find that the ratio f a dv-g/ fheat-g exceeding 
unity corresponds roughly to the time when the shock begins 
its expansion towards explosion in two-dimensional runs. Ob- 
viously all the saturation mechanisms discussed in jj4.5l are 
present in those simulations, so analogi es need to be made 
cautio usly. A similar result is obtained by Murphy & Burrowsl 
(2008), with a heating time defined without gravity and ki- 



netic energy. 

An attempt to include the effects of convection in the 
steady -state solution was made by Yamasaki & Yamada 
(2006), setting the convective flux to a value that yielded a 
flat entropy gradient. With this maximally efficient convec- 
tion, they found that the Burrows-Goshy limit can decrease 
by several tens of percent relative to the laminar case. This 
would entail a corresponding decrease in the approximate 
instability thresholds found in this paper. A more careful 
treatment of the turbulent transport terms, along the lines of 
Murphv & Meakin (2011), could be used to construct quasi- 
steady-state flows and to per form a linear analysis like that of 
Yama saki & Yam ada (2007), helping to elucidate these ques- 
tions with semianalytic tools. 

The extension of our model to two- and three spatial dimen- 
sions will be discussed in a companion paper. 
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APPENDIX 
A. WEAK INTERACTION RATES 

Here we describe the calculation of weak interaction source terms entering the energy and lepton conservation equations 
(eqns. H and Q, respectively). We only consider electron-type neutrinos and antineutrinos, as other species do n ot exchange 
energy with matter outside the neutrinosphere. Our implementation largely resembles that of Ohnis hi et al.l ([2006), with some 
minor modifications. 

For both species, we assume a free-streaming neutrino flux radiated isotropically from neutrinospheres located at the same 
radius R u . The neutrino distributions are approximated by a Fermi-Dirac function with zero chemical potential 

f Vi (e, T Vi , r, cos k ) = N Vi F FD (e, T Vi , 0) 9 [cos k - cos 6 v (r)] , (Al) 

where v i = {v e ,v e \, 

Ffd(^ T, /i) = \ (A2) 

exp(e-/i)/fer+ 1 

is the Fermi-Dirac distribution with energy e and chemical potential /i. The normalization factor 

= (7/16)47^^4 (A3) 

is such that the neutrino luminosity L v , and spectral temperature T Vl are independent parameters. The step function 9 contains 
the angular dependence of the radiation field, which propagates up to an angle 

1/2 

(A4) 



0k < Qi>(f) = arccos 



from the radial direction. In equation ( IA3I ). ctsb is the Stefan-Boltzmann constant. 
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To approximate the results of radiation-hydrodynamic simulations with our simplified assumptions, we set L Vc = L Pe but specify 
different temperatures, T Ve = 4 MeV and Tp e = 6 MeV, reflecting th e fact that the pho tospheres are at slightly different locations, 
but with a spectrum such that luminosities are roug hly comparable dJankall 19951 [2001 ). We ignore however the correction due to 
the spatial difference between the electron neutrino- and antineutrinospheric radius, which is of the order of 10%. This difference 
is absorbed by the normalization factor in equation ((A3). 

The evolution of Y e above the neutrinosphere is determined by the net production rate of electron and positrons IV and r e +, 
respectively (eq. JTJ), due to neutrino absorption and emission. Following Bruenn ( 1985), these are given by 



2Ttm„C 

r, 



(he) 3 
27rm n c 



(hc) 3 p 



p J dcostf* J e*de [K u J Ue -j Ve (l-fu.)] 
2 de , 



dcos#t 



(A5) 
(A6) 

(A7) 



where j Vj and K Vi are the emissivity and absorption coefficient, respectively, associated with electron-type neutrinos or antineutri- 
nos (as subscripted), and the integrals ar e performed over all propagation angles and positive energies. Expressions for these co- 
efficients are obtained by Bruenn ( 1985) assuming detailed balance, matter in nuclear statistical equilibrium, and non-relativistic 
nucleons, whose recoil is neglected. In addition, we ignore here the nucleon phase space blocking factors, which deviate only 
slightly from unity at densities < 10 11 g cm" 3 dBruennll 1985b . yielding 



ju,(e,T, fi e ,n p )-- 
jff e (e,T,fj, e) n p )-- 



(g 2 v + 3g 2 A ) n p F FD (e,T,p e )[e+A m f x 



- (gv+3gl) n„F FD (e-A m ,T-n e )[e-A m ] 2 



2 4 



(e + A m ) 2 



1/2 



7T 

6(e-A 



2 4 1 V2 



(e-A m ) 2 



-m e c 2 ) 



K Ve (e, T, fi e , n p ) = — (gy + 3g 2 A ) n n [1 - F FD (e , T, fi e )] [e+ A,„] 2 



2 4 



(e + A,„) 2 



1/2 



G 2 

Kv e (e, T, (i e ,n p ) = — (g 2 J + 3g 2 A ) n p [1 -F m (e- A„„ T,-fj, e )] [e- A m f 



7T 

6(e-A 



2 4 



(e-A m ) 2 



1/2 



-m e c 2 ), 



(A8) 

(A9) 
(A10) 

(All) 



where /j, e is the chemical potential of electrons, n„ and n p the number density of free neutrons and protons, respectively, G F = 
Gf /(he) 3 the Fermi constant, gy and g A the vector and axial coupling constants, respectively, and A m = {m„-m p )c 2 the difference 
between the rest mass energy of neutrons and protons. 

The rates of energy exchange between neutrinos and matter are similarly found as a function of f Vi , j Vi , and n Vi . The heating 
and cooling rates per unit volume are given respectively by 



(A12) 
(A13) 



where are rates per baryon of heating (+ superscript) and cooling (- superscript) due to electron-type neutrinos {y e subscript) 
and antineutrinos (p e subscript). These rates are given by 



2irm n c 
(hc) 3 p 
2nm n c 

' (hc) 3 p 
Anm n c 

' (hc) 3 p 
4irm n c 

: (hc) 3 p 



d cos 9 k J e 3 d e [j Ue + k„J /„, 
d cos 9 k J e 3 d e [j Pe + K Pe ] f Pe 



(A14) 
(A15) 
(A16) 
(A17) 



Numerical calculation of equations ( IA6t -( tA~7b and (IA14b - (lAT7l involves tabulating Fermi-Dirac integrals as a function of T, 
H e , and T Vj , with other dependencies entering as global scaling factors. The left panel of Figure [A] shows the total cooling per 
baryon Q^ ot = Q~ +Qp as a function of temperature, for fixed Y e and different densities. In the degenerate regime kT < irp e , the 
cooling rates are nearly independent of temperature, whereas for higher temperatures they approach the asymptotic expression 
Q~ ol ~ 145(A;r/2MeV) 6 from ljaiikal(r200Tl) . obtained by assuming zero electron chemical potential. The right panel of Figure [A] 
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FIG. 16. — Left: Total cooling rate per baryon Q[ oi = Q e _+Q e+ (eqns. IA16I - IA17I ) as a function of temperature, for fixed electron fraction Y e = 0.3, and assuming 
X p +X„ = 1. Different solid lines correspond to different densities (in g cm -3 ), with circles marking the approximate boundary between degenerate and non- 
degenerat e regimes, at /i e = irkT. The red dashed line shows the approximate expression 0^, — 145(W/2MeV) 6 MeV s , obtained assuming non-degenerate 
electrons I Janka 2001). Right: Heating due to absorption of electron-type neutrinos (eq. IA14I ) as a function of the temperature of the electron neutrinosphere 



T Vr , for para meters as shown in the panel. For reference, the red-line shows the approximate expression Q^ e 
IjanklEOOlb . 



160^,52^(1- 



Y e )(kT v J4MeVf MeV s" 1 from 



shows heating by absorption of electron-type neutrinos as a function of T Ve , for fixed parameters as shown in the Figure. For 
reference, the approximate expression Ql e ~ 160L^^2r^ 2 (l-Y e )(kT^j4-MeW) 2 MeV s" 1 (Janka 2001) is also shown, agreeing 
with our results to within ~ 10%. 

B. GRID OF VARIABLE SPACING IN FLASH3.2 

Time-dependent modelling of post-bounce core-collapse supernova hydrodynamics demands the ability to resolve a steep 
density gradient in the protoneutron star atmosphere (~ 30 km) while also allowing for the shock to expand out to at least 
~ 1000 km to track an explosion. For a grid in spherical coordinates, an efficient way to accomplish this is to use variable 
spacing in radius. 

The public version of FLASH3.2 was modified to include this capability, starting from the existing uniform grid mode. The 
implementation can be decomposed into two parts. First, defining the cell sizes and coordinates appropriately when the compu- 
tational domain is initialized. The second part involves modifying all the subroutines that assume a uniform grid spacing. 

Grid initialization is accomplished by modifying the Grid_init and gr_create_domain subroutines. We define the 
grid points in between r m [ n and r max such that consecutive cell sizes have a ratio Ar l+ i / Ar; = C, > 1, where i is a cell index which 
increa ses with increasing radius. For a given number of grid cells N r the ratio can be obtained by solving (e.g. JStone & Norman! 
[19921) 



P max ^ m j n — miQ ^ ^ C 



■■Ar 



C N >-i 



c-i 



(Bl) 



where Ar m ; n is the minimum cell size, whose inner edge is located at r m i n . A logarithmic spacing can be achieved by setting 
C = (''max /'"min )* , which then determines the minimum cell size: Ar^ = r mm (£- 1). It then follows that for < q < N r , 



C" 



Ar n 



Ar n 



iXql ' min) 



(B2) 



with r = r min and r Nr = r max . 

There are only three subroutines that make the assumption of a uniform grid when using spherical coordinates. They are 
hy_ppm_sweep, Driver_computeDt, and Driver_verif ylnitDt. In all three cases, a scalar cell spacing was re- 
placed with a vector in the appropriate locations. The subroutines that compute cell areas and volumes make direct use of 
coordinate information, so no additional modification is required for spherical coordinat es. 

As a test of the grid implementation, we have run the self-similar explosion problem of lSedovl ([1982). An energy Eq is initially 
placed inside some spherical volume with a radius smaller than some characteristic radius Rq. The medium has uniform density 
po, and has an ideal gas equation of state with adiabatic index 7. Outside of the injection volume, the pressure is uniform an equal 
to lO^Zio^o 3 - O ur benchmark simulation has 500 cells uniformly spaced in radius from r = to r = Rq. The energy is placed in 
the first grid cell, with radius r; n j = 2 x 10~ 3 /?o, and we set 7 = 1 .4. 
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FIG. 17. — Strong explosion problem of Sedov 1 1982), used here as a test of the non-uniform grid implementation in FLASH3.2. The left panel shows the 
shock trajectory as a function of time, with colors labeling resolutions as shown in the caption. The minimum cell size next to the origin is kept equal to the 

-1/2 1/2 5/2 

uniform grid case = 1). The right panel shows the fractional deviation from the shock position at time 0.75£p p R in the N r = 500 ran, with error bars 
indicating the cell sizes in each simulation at their respective shock positions, normalized to the shock position in the N r = 500 run. 



We then run a sequence of simulations that preserve the minimum cell spacing next to the origin, Ar^n = r^, but vary 
the number of cells, resulting in a ratioed (not logarithmic) grid that satisfies equation ( IBU . The number of cells are 
N r = {466,400,300, 180} which results in (£- 1) ~ {3 x 10" 4 , 10" 3 ,3 x 10" 3 , 10" 2 }, respectively. The left panel of Figure [B] 
shows the shock radius (obtained by linear interpolating the location of the surface with pressure 10~ 3 EqR^ 1 ) as a function of time 
for the sequence of runs. As the number of cells is decreased relative to the benchmark model, the cells in the upper part of the 
domain become increasingly coarser, and the shock trajectory gradually diverges from the uniform grid case. The right panel of 

FigurelBlshows the fractional deviation from the uniform grid result at a time 0.75Eq p Q R , with the error bars indicating the 
fractional size of the cell at the given shock position, normalized to the shock position in the benchmark model. As the spacing 
ratio decreases, deviations tend to zero close to linearly in (£— 1). 

C. NUMERICAL CALCULATION OF THE RATE OF CHANGE OF ENERGY 

Here we describe the calculation of the different terms that comprise the rate of change of the total energy with time (eq. fl9l ) 
for Figures [5] and [7] 

In the Piecewise Parabolic Method (IColella & Woodwardll 19841) used in FLASH, the shock is typically broadened along 2-3 
cells. To compute an accurate shock position, we begin by finding the minimum in the radial derivative of the pressure. We 
establish a reference position by taking a weighted average of the radial coordinate, with weight equal to \dp/dr\ around the 
minimum. This reference position, which varies smoothly with time, is a good tracer of the center of the shock. Because we want 
quantities below the shock, we define our actual shock position to be a fixed distance (of the order of two cell widths) behind 
the center of the shock, so that the fluid quantities are in the post-shock regime while remaining as close to the shock center as 
possible. 

Next, the energy flux entering through the shock E up (eq. ETlP is computed by interpolating variables (linearly for velocity and 
internal energy; logarithmically for pressure and density) to the shock position. The energy flux leaving the domain E^ n (eq. 1221 ) 
is computed at a radius Ri n corresponding to the inner edge of some cell inside the domain. After some experimentation, we found 
that the innermost location that is not significantly affected by boundary effects is the third active cell from the inner boundary at 
R„. To compute the energy flux at we average the fluid variables from the cells that this radius separates as an approximation 
to face-centered values. The integrated neutrino source terms (eq. j20l ) are computed by simple integration of the cells fully 
contained between R\ n and the shock, and then a small differential term is added by linear interpolation, 47r/? 2 mt J;? net AR, where 
i? ut is the outer radius of the cell that is closest to and fully contained within the shock surface, and AR = Rs~R ut- The shock 
motion term E$ is found by first computing the shock velocity through time-centered finite differences, and then evaluating fluid 
quantities at the same position as for E up . 

For calibration, we also compute the total energy contained between R m and the shock. The total rate of change, computed 
through time-centered finite differences, and the sum of the separate terms that make it up (eq. fl9l ) are shown in Figure [18] 
Quantities are sampled every 2 ms to eliminate high frequency noise. Agreement is very good if a zero point, equal to the sum of 
all the terms that make up E tot at t = 0, is subtracted. 

Also shown in Figure[l8]is the net energy generation from accretion, E^-E up +Ei n , together with what is obtained if the shock 
motion term E s is subtracted from the total rate of change of the energy E tot , computed from finite differencing the total energy. 
Agreement is again very good if the same zero point is subtracted. This smoother version of the accretion energy generation, 



HYDRODYNAMICS OF CORE-COLLAPSE SUPERNOVAE 



21 



10 



S> 5 



s o 



■a 

+ 

■uf 

' -5 

•w 



_ I I 1 , 1 , , , , 1 1 , 1 , 1 , , , 1 1 1 1 1 1 1 

— sum of terms 


i i i i 1 i _ 


~ — finite difference 








i i i i i i i i i i i i i i i i i i i i 


III 

/X" 


U — ^-^vaA 




. ... 1 . ... 1 .... 1 .... 1 .... T 


.... 1 . 







0.05 



0.1 



0.15 



0.2 



0.25 



0.3 



time [ms] 

FIG. 1 8 . — Top: evolution of the total energy in the model shown in Figure[5] The black curve corresponds to the sum of the terms on the right hand side of 
equation )19t . calculated as described in the text. The red curve is a time-centered finite difference derivative of the instantaneous total energy. Bottom: evolution 
of the net energy generation from accretion. The black curve is the sum of the terms that comprise it, while the red curve is obtained by subtracting the energy 
generation from shock motion E s from the red curve in the top panel. 

together with E to t from finite differencing, is what is shown in Figures [5] and [7] 
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